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Abstract. Dynamically cold components are well known to destabilize hotter, even much more massive components. E.g. 
stellar disks can become unstable by a small admixture of cold gas or proto-planetary disks might be destabilized by a small 
fraction of dust. In this paper we studied the dynamical influence of a cold dust component on the gaseous phase in the central 
regions of galactic disks. We performed two-dimensional hydrodynamical simulations for flat multi-component disks embedded 
in a combined static stellar and dark matter potential. The pressure-free dust component is coupled to the gas by a drag force 
depending on their velocity difference. 

It turned out that the most unstable regions are those with either a low or near to minimum Toomre parameter or with rigid 
rotation, i.e. the central area. In that regions the dust-free disks become most unstable for high azimuthal modes (m ~ 8), 
whereas in dusty disks all modes have a similar amplitude resulting in a patchy appearance. The structures in the dust have a 
larger contrast between arm and inter-arm regions than those of the gas. The dust peaks are frequently correlated with peaks 
of the gas distribution, but they do not necessarily coincide with them. Therefore, a large scatter in the dust-to-gas ratios is 
expected. The appearance of the dust is more cellular (i.e. sometimes connecting different spiral features), whereas the gas is 
organized in a multi-armed spiral structure. 

We found that an admixture of 2% dust (relative to the mass of the gas) destabilizes gaseous disks substantially, whereas 
' dust-to-gas ratios below 1% have no influence on the evolution of the gaseous disk. For a high dust-to-gas ratio of 10% the 

instabilities reach a saturation level already after 30 Myr. The stability of the gaseous disks also strongly depends on their 
q . Toomre parameter. But even in hot gaseous disks a destabilizing influence of the dust component has been found. 
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1. Introduction Revealing the nature of the mini-spirals in the central re- 

■ gions of galaxies is of great importance for our understanding 

Recent high resolution optical and sub-mm observations of the of a variety of astrop h y sical processes such as the mass ac- 

central regions of galaxies have revealed that the circumnuclear cretion rate of the central en gines of galaxies or the evolution 

regions contain stellar-gaseous mini-disks with a size of about of the ga i actic disks themselves. Therefore, the dynamics of 

a few hundred pc An extensive observational program on HST gase0 us density waves in the central region of disk galaxies 

by Carollo et al. Umi\mM\im^ who studied the nuclear has attracte d more a nd more attention during the last decade, 

regions of seventy five early and intermediate type disk galax- Athanassoula i fT^ studied the gas flows in and around bars, 

ies, unveiled an astonishing richness of structure in the nuclear she demonstrated by orbital analysis and 2D-hydrodynamical 

regions of galaxies. About sixty percent of their sample have simu i ations th at the dynamical properties of the disk in connec- 

mini-bars, spiral-like dust lanes, star-forming rings and spiral tion with a rotating large-scale bar, especially the existence of 

arms. Many nuclear regions are also well described as patchy different periodic orbit families (or equivalently the existence 

or mul ti-armed. Further investigations by Regan & Mulchaey of inner Lind blad resonances (ILR)) are key criteria for the ex- 

( 1999 1 showed that spirals are the most common morphologi- istence and shape of dust lanes If X2 and ^-orbits are exist- 

cal structures in the central regions of galaxies. Often the spiral ing> gas (and dust) follow first X[ . orbitSi until they enter 

patterns are not associated with the outer grand-design spiral regions which support also " anti _ bar " orbits . In these turbulent 

patterns of these galaxies which points to different physical ori- reg i ons the gas dissipates energy and switches to x 2 or * 3 -orbits 

S lns - which brings the gas to the central area. 

Different to the interpretation based on a stellar orbital anal- 
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spirals in the central regions of galactic disks are related to the 
formation of grand-design spiral patterns in galaxies. They ar- 
gued that gas density waves - different to stellar density waves 
- are not completely damped or absorbed at the ILR and, thus, 
they may generate spiral structures at all radii including the nu- 
clear regions. Such a model might explain the continuity of the 
spiral features at small and large radii as well as the low arm- 
interarm contrast observed in galactic centers. This scenario is 
supported by their numerical simulations studying the gas re- 
sponse to a background gravitational potential composed of a 
galactic disk and a large-scale stellar bar. 

An alternative application of barred potentials as driving 
engines for the generation of nuclear spirals is to invoke sec- 
ondary bars or small nuclear bars located well inside the ILR. 
Wada & Koda (2001) performed a set of multi -phase hydro- 
dynamical simulations studying the dynamics of central mini- 
disks influenced by the potential of a weak mini-bar. They 
took into account the self-gravity of the gas as well as cool- 
ing and heating processes. They found the formation of fil- 
aments and cusps on a parsec-scale which resemble the ob- 
served morphological patterns in the central regions of galax- 
ies. Another good example for the influence of a nuclear bar is 
the match between observations and numerical (SPH-)models 
by Ann (2001 1. He reproduced well the shape and orientation 
of the nuclear ring of the barred galaxy NGC 4314. Jogee et 
al. (|2002 1 studied the propagation of density waves triggered 
by bar shocks which were derived from multiwavelength ob- 
servations of NGC 5248. Motivated by the existence of mas- 
sive molecular arms inside the disk they incorporated the self- 
gravity of the gas. Their comparative analysis shows an agree- 
ment between the modelled and observed gas morphology, gas 
kinematics, and pitch angle of the spirals. These results argue 
for a dynamical coupling between the nuclear region and the 
surrounding disk. 

A common property of all the mentioned mechanisms is 
that they usually result in structures dominated by two or a few 
arms. E.g. a nuclear bar would naturally induce a two-armed 
structure. Similarly, grand-design features penetrating from the 
outer disk through the ILR will also keep their symmetry prop- 
erties characterized by low azimuthal mode numbers. However, 
Elmegreen et al. (|2002 1 found that nuclear dust spirals differ 
from main-disk spirals in several respects: the nuclear spirals 
are not associated with star formation, they are very irregular 
with both trailing and leading components that often cross, and 
they completely fill the inner disk with a constant surface den- 
sity. 

Little is yet known about the dynamics on the 10-100 pc 
scale of the central regions of galaxies. Attempts to explain the 
existence of such structures as a result of gravitational instabil- 
ity of the dynamically decoupled mini-disks were believed to 
be unsuccessful because of the low surface densities of the fast 
spinning gaseous mini-disks. As an alternative to gravitational 
instability as an origin of structure in mini-disks Elmegreen 
et al. (1998 1 suggested that mini-spirals may be the result of 
an acoustic amplification of the grand-design spirals propagat- 
ing to the central regions of galaxies. Montenegro et al. ([1999 1 
demonstrated that such an amplification mechanism can indeed 
work inside the ILR or outside the OLR. 



All discussed mechanisms have some weak points and can 
not fully account for all the observational data. We consider a 
new approach to understand the origin of mini-spirals in the 
central regions of galaxies. Observations show that the circum- 
nuclear disks of galaxies are dusty, and dust may play an impor- 
tant role in their dynamics. It is known that cold components 
greatly destabilize multi-component gravitating disks (Jog & 
Solomon 1984] Bertin & Romeo[l988 Orlova et al. l20U2l . An 
admixture of a dynamically cold dust component to the gaseous 
mini-disks of galaxies might therefore strongly destabilize the 
disk leading to the formation of spiral structure. From a stabil- 
ity analysis Noh et al. ( 1991 1 showed that a dust component can 
strongly destabilize proto-planetary disks. The admixture of 
only 2% of dust can enhance the growth rates of the dominant 
gaseous phase significantly. This effect becomes even larger 
with less massive disks. A conservative estimate of the dust- 
to-gas ratio in the solar neighborhood gives an average value 
of 0.6% ranging from 0.2% up to 4% in H2 regions (Spitzer 
I1978> . These are lower limits because large grains which con- 
tain much mass do not have any detectable extinction. Recent 
observational data by Maiolino et al. (2001 1 exhibit evidence 
for anomalous properties of the dust grains in the galactic nu- 
clei. By comparing the reddening of optical and infrared broad 
lines and the X-ray absorbing column density they found that 
the Ay/Nn ratio is significantly lower in the circumnuclear re- 
gions of galaxies than in the diffuse regions of galaxies. This in- 
dicates that the dust composition in the circumnuclear regions 
of galaxies could be dominated by large grains and, thus, sub- 
stantial amounts of dust might have been undetected so far. 

In this paper we study whether the origin of mini-spirals 
can be explained by a destabilizing role of a dust component 
in the circumnuclear disks. We restrict our analysis here to sys- 
tems without a nuclear or a large scale bar. Also we do not 
consider here any formation or destruction processes of dust. 
Thus, we investigate the impact of a frictional force exerted 
by the interaction of dust and gas on the dynamics of galac- 
tic nuclear regions. This investigation is done by means of 2- 
dimensional multi-component hydrodynamical simulations. In 
the next Section we describe the numerical model, i.e. the basic 
hydrodynamical equations and the dust-gas interaction as well 
as our numerical code. In Sect. [3] the results of the numerical 
simulations are shown which are discussed in more detail in 
Sect.@] Finally, a summary is given in Sect.|5] 

2. Numerical method 

2.1. Pure hydrodynamics 

We studied numerically the hydrodynamical equations for a 2- 
dimensional single- or multi-component disk. Thus, we solved 
the continuity equation 

and the momentum equations for gas and dust. The momentum 
equations read for gas 

^ + (v s • V)v g + ^ + V(0> + Ohbsd) = S g (v g ) (2) 
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and for dust 

^f+(Vd- V)v rf + V((D + 0>hbsd) = S d (v d ) ■ (3) 
at 

The components are denoted by g for the gaseous phase and 
by d for the dust. Y, gtd are the surface densities and v gtd the 
velocities. P ? is the pressure of the gas which is given in the 
2d-case as force per unit length. cD denotes the potential of the 
self-gravitating disk. It is derived from the Poisson equation 

AO = 4nGI.(R, <p) 6(z) = AnG (s g + S d ) 6(z) . (4) 

An external potential is added by a stationary contribution 
Ohbsd related to the halo, the bulge and/or a stellar disk compo- 
nent. These external potentials are chosen to match - together 
with the potential of the disk - a given rotation curve. 

The main difference between the gaseous and the dust 
component is that the dust is treated as a pressureless phase. 
Therefore, if gas and dust are in rotational equilibrium, there 
is a velocity difference between both components which might 
give rise to a non-negligible frictional force depending on the 
cross-section for the gas-dust interaction. This interaction is de- 
scribed by the source terms S(, , .) on the RHS of the hydrody- 
namical equations. Since we do not consider dust formation 
and destruction processes, the source terms in the continuity 
equation vanish. However, frictional terms show up in the equa- 
tions of motion. The dust implementation will be described in 
the next paragraph in detail. 

The set of hydrodynamical equations is closed by a poly- 
tropic equation of state 



KZ? g s 



(5) 



For the gaseous phase we set the polytropic exponent to y g - 
5/3. The constant K is chosen to yield a given minimum 
Toomre parameter. 

2.2. Treatment of the dust component 

Gravitating mini-disks in galaxies differ from "normal" galac- 
tic disks in which the cold component (gas) is coupled gravi- 
tationally to the dynamically hot component (stars). The main 
difference between galactic disks and mini-spirals is that the 
cold component (dust) in the mini-disks is not only coupled by 
gravity to the hotter component (gas), but also by a frictional 
force between both components. 

In order to include this drag, we added a source term to the 
equations of motion following the general form suggested by 
Noh et al. 



/ = S d (v d ) 
S g (v g ) 



-A(v d - v ? ) 
2. 



(6) 
(7) 



The second source term, Eq. Q follows from the requirement 
of momentum conservation. The physics of the friction is en- 
closed in the frictional timescale A -1 . We implemented three 
different prescriptions of the frictional term in our models. 

The simplest approach is to assume a time- and position- 
independent frictional timescale r d , i.e. 
-l 



Though this is a very rough assumption, it allows for a direct 
comparison between the drag term and other dynamical quanti- 
ties, e.g. the dynamical timescale given by the rotation period. 
More detailed approaches are based on a microscopic view of 
the frictional process, i.e. the momentum exchange between 
gas and dust particles by collisions and the equipartition of mo- 
mentum within the gaseous phase. Two limiting situations are 
described in the following sections in more detail. 

2.2.1 . Collisional time scale (equal disk heights) 

Because the size R d ~ 0.1/mi of typical dust particles is much 
smaller than the mean molecular free path A ~ 0.1 AU of the 
ambient gas particles in galactic central regions, the veloci- 
ties of gas molecules are uncorrelated with the velocity of the 
dust particle itself and, thus, the friction can be described as an 
"Epstein" drag (Goodman & Pindor 2000 1. The latter is derived 
from kinetic gas theory (Epstein 119241 . Only for very large 
grains residing in very dense regions the Stokes' formula for 
the friction has to be used (then the friction is given by a lam- 
inar viscous flow over the dust particle). Assuming "Epstein" 
friction is acting Noh et al. ( 1991 1 discerned between the two 
limiting cases of different scale heights of the dust component: 
either similar disk heights of gaseous and dust disk or a much 
thinner dust disk. 

In case of similar or equal disk heights it is a reasonable 
ansatz to use the time scale t c for energy and momentum ex- 
change due to collisions between gas and dust particles as an 
estimate of A -1 
. _i o- c p g v th 



m d 



(9) 



md denotes the mass of a dust particle, p g the mass density of 
the gas and v t h the thermal velocity of the gas. The collisional 
cross section cr c can be estimated from the geometrical cross 
section, i.e. from the average radius R d of a dust particle 



cr c =n-R d . 

The collisional time scale is then given by 
r, * 1.49- 10 3 



(10) 



m d \{ Rd 
10- 14 g/\0.1/mi 



H 



10 3 M©pc- 2 / \100pc/ \l0kms 



Vth 



yrs (11) 



H is the full scale height of the gas, its spatial density is esti- 
mated by p g = "LglH. 

2.2.2. Dynamical timescale (thin dust disk) 

The other limiting case corresponds to a dust disk which is sub- 
stantially thinner than the gaseous disk. In that case the time 
scale should be of the order of the dynamical time scale given 
by the inverse circular frequency Q. 1 (Noh et al. U991l l. i.e. 



A - Q. 



(12) 



a = t; 



(8) 



The basic assumption is here that gas and dust establish very 
fast collisional equilibrium in the thin dust layer and that the 
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gas momentum is then mixed vertically in a sound travelling 
time scale. When taking typical values we see from Eq. il Q 
that the collisional time scale t c is indeed much shorter than 
the dynamical timescale. E.g. the rotation period in the inner 
regions of galaxies is about 10 7 yrs. The sound propagation 
time t s can be estimated from t s ~ H g /v^ (H g is here the disk 
scale height of the gaseous phase). In dynamical equilibrium 
the scale height is related to the velocity dispersion by H g ~ 
v 2 h /(jrG"E) (see e.g. Binney & Tremaine 1987 hereafter BT87) 
which yields 



Vth 

ttGS 



cr 



The second relation is derived from the Toomre stability crite- 
rion. The last relation between the circular frequency Q. and the 
epicyclic frequency k holds for all reasonable rotation curves 
within a factor of 2. Hence, in the limit of a thin dust disk the 
friction time scale is dominated by the sound propagation time 
which is of the order of the dynamical time. 

2.3. Numerical implementation 

The nonlinear analysis implies the solution of the full set of hy- 
drodynamical equations. For this purpose, we developed a two- 
dimensional numerical code which is similar to the ZEUS-2D 



code by Stone & Norman ( 1992 1. The hydrodynamical equa- 
tions are discretized on a Eulerian grid in polar coordinates. 
The different terms are treated by operator splitting. Advection 
is performed by a second order Van Leer advection scheme. 

In radial direction we use a logarithmic grid which allows 
a very high resolution in the central region. Therefore, a ra- 
dial resolution of 128 cells should be sufficient for two-armed 
spirals as Englmaier & Shlosman (2000 1 demonstrated. In az- 
imuthal direction, however, 128 grid cells can be critical, be- 
cause the maximum growth rate might be found for high val- 
ues m of the azimuthal wavenumber. In case of protoplanetary 
disks Noh et al. {1991 1 found m ~ 5 ... 10 to be the most un- 
stable modes. In order to accommodate such high modes, we 
used mainly a grid size of 270x270 cells 1 . The physical sizes 
of the cells are similar in azimuthal and radial directions at all 
distances. The typical radial extent of the whole grid covers 
radial ranges from a few 10 pc to a few kpc. The Poisson equa- 
tion is solved by applying the two-dimensional Fourier convo- 
lution theorem in polar coordinates (BT87). The hydrodynam- 
ical timestep is calculated by a combination of the standard 
Courant-Friedrichs-Levy criterion and the timescale is derived 
from the frictional force. 

2.4. Units 

For our simulations we set the gravitational constant to unity 
and we chose the length unit to 1 kpc and the mass unit to 

1 For our simulations on the NEC-SX5 of the Kiel computing cen- 
ter we used a 270x270 grid which has a similar physical resolution 
as the "normal" 256x256 grid, but a much better computational per- 
formance of the vector-optimized FFT routine supplied with NEC's 
MATHKAISAN library. The increase of performance is caused by 
avoiding bank conflicts in the memory access on the NEC-SX5. 



10 9 Mq. The system time is then measured in 1 .49 • 10 7 yrs, the 
velocity in 65.6 kms~', the circular speed in 65.6 kms -1 kpc -1 
and the surface density in 10 3 Mq pc -2 . If not stated explicitly, 
all quantities are given in these units. 



2.5. Initial models 

For most of our simulations the mass distribution of the disk is 
described by a weakly perturbed exponential profile. The per- 
turbation is created by multiplying the surface density with a 
factor (1 + Cr). r is a random number in the interval [-1,1] and 
the amplitude C is typically selected to yield a global Fourier 
amplitude of the order 10~ 6 (cf. also Eq. J18l l). 

The initial velocities are derived from rotational equilib- 
rium, i.e. there is no initial radial motion. The rotation curve is 
specified by 

R 

■Rflat 



V C {R) - V'o. 



1 + 



R 



l/n, 



(13) 



The radius Rn- dt defines the transition between the central rigid 
rotation and an outer flat rotation curve. The sharpness of the 
transition is controlled by n,. The azimuthal velocity of the 
gaseous phase is calculated by the (frictionless) Jeans' equa- 
tion for the radial velocity component which reads in cylindri- 
cal coordinates 



du g 

+ w „ 
dt dR 



du g v g du g 



R d(f> 



_£ 
R 



1 



dR 



d 

— (O + <3>hbsd) 
dR 



.(14) 



In case of the initial equilibrium, the radial velocity m ? vanishes 
and, hence, the azimuthal velocity v g is given by 



v g 2 = 



R dPg 
ET ' ~dR 



+ R— (O + Ohbsd) = v 2 c (R)-v 2 P (R) . (15) 



where we used that the rotation curve is derived from the gravi- 
tational potential, i.e. v 2 (R) = Rjfft(Q>+Q>HBSD), and the pressure 
contribution Vp to the rotation curve is abbreviated as 



r_ 

S e ' dR 



(16) 



In case of the pressureless dust component its initial azimuthal 
speed is directly given by the rotation curve dl3l >. 

From the mass and potential distribution (or the rotation 
curve), the radial shape of the profile of the Toomre parameter 

Vth* 



(17) 



can be derived except for a constant factor. This factor can be 
specified by requiring a minimum value Qmin f° r the Toomre 
parameter. 

2.6. Fourier Modes 

In order to quantify the (global) stability of the disk we use the 
global Fourier amplitudes of each component (cf. e.g. Laughlin 
et al. lT998l 

1 



r = 

— 



M, 



disk 



rvLll r-R^ 
JO jR m 



2(r, <p)rdr e 



(m > 0). 



(18) 
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Fig. 1. Temporal evolution of the Fourier amplitudes of the m 
0, 1,2,6,8, 10-modes for a single component model with the 
mass distribution of the gaseous component in the reference 
model. The time unit is 1.5 ■ 10 7 yrs. 



Mdi s k is the mass of the disk for the component of interest in the 
specified radial interval [R\ n , R out ]. S denotes the corresponding 
surface density. Global Fourier amplitudes are calculated by 
integrating over the whole disk. 

Additionally, a mode m — is introduced which quantifies 
the radial mass redistribution with respect to the initial one: 



C = 



In 
Mdisk 



L 1 



\t(R, t) - 2o(R)\RdR 



(19) 



£(/?, t) and %)(R) are the azimuthally averaged surface densities 
at time t and at the begin of the simulations, respectively. As 
long as the evolution is well described by linear perturbation 
theory, no substantial radial mass transport occurs. Thus, the 
increase of Co is a good indicator for the onset of non-linear 
effects during the growth of perturbations. 

2. 7. A single component model 

In order to study the properties of a standard disk, but also 
to test the code for our application we performed a series of 
single-component (i.e. purely gaseous) simulations based on 
the mass distribution of a reference model. Its mass distribu- 
tion is characterized by a submaximal exponential disk and a 
rotation curve which increases linearly in the center and is flat 
outside (for exact values see the next section). 

In order to quantify the structural evolution we first dis- 
cuss the global Fourier amplitudes (Fig. [3: they remain con- 
stant on their initial level of about 1CT 6 during the first 100 
Myrs (or t » 7). Then the high-order modes m = 6,8,10 
start to grow exponentially reaching a saturation level of about 
10~ 3 . The fastest growing and dominant mode is the m=8- 
mode. This is also seen in the image of the density perturba- 
tions (Fig. |2 upper row). The wavelengths of the perturbations 
are rather short as expected from the small critical wavelength 
A c = 4n 2 GL/K 2 x 100 pc. From linear theory the maximum 
growth is expected for a wavelength at about A c /2 ~ 50 pc 
(see BT 87). Later than the high order modes, the amplitudes 
of the low-order modes, especially m = 2 and m = 1, com- 




mamm. 



1& 




Fig. 3. Temporal evolution of the spatially resolved Fourier am- 
plitudes for the ;7i=2-mode (upper diagram) and the m=8-mode 
(lower diagram) for a single component model with the mass 
distribution of the gaseous component in the reference model. 
The amplitudes are calculated according to Eq. Jl 8i integrating 
over small radial annuli defined by the numerical grid resolu- 
tion. The time unit is 1.5 ■ 10 7 yrs. 



mence to increase with a time delay of 100 Myr. At the end 
of the simulation at t = 40 (600 Myr) their saturation levels 
slightly exceed the saturation levels of the high order modes. 
The overall appearance is then dominated by a mainly lopsided 
and very irregular morphology (Fig. [2 lower row). 

The radial mass distribution varies during the first 600 Myr 
only within the innermost 100 pc which is reflected also in the 
evolution of the m=0-mode: After a fast readjustment of the ra- 
dial mass profile according to the initial density perturbations, 
it remains constant until t ~ 15 » 225 Myr. Then a radial 
mass transfer sets in, which stops at about t ~ 30 w 450 Myr. 
Though the values of the perturbation amplitudes are rather 
small, this does not mean that the response of the disk is ev- 
erywhere in the linear regime. The outer regions remain almost 
unevolved, whereas the inner regions undergo strong morpho- 
logical changes. It is only the "standard" normalization to the 
whole disk mass which keeps the global Fourier amplitudes 
rather small. In order to demonstrate this, Fig.[5]shows the tem- 
poral evolution of the spatially resolved Fourier amplitudes for 
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Fig. 2. Spatial distribution of the surface density perturbations (normalized to the their initial values) of the single component 
gaseous nuclear disk (Qmin = 1.54) at different times: during the linear growth regime (f = 16, upper left), saturation of the 
higher-order moments (f = 24, upper right), saturation of the lower-order moments (f = 32, lower left) and at the end of the 
simulation (t = 40 ~ 600 Myr, lower right). Areas devoid of material are white, areas with a density enhancement of a factor of 2 
or more are black. The grey area seen at the outer edges correspond to no deviation from the initial surface density. Note that the 
data inside the inner boundary of the computational grid (|jc|, |y| < 0.03 = 30 pc) are artificial. They are created by the graphics 
software when converting from the polar grid to the Cartesian grid of the image. 



m-2 and m—%. As seen already for the global modes, the high- 
order modes start to grow first. Additionally, the growth be- 
gins in the central area where the surface densities are at max- 
imum and the dynamical timescales are short. Throughout the 
whole simulation the region of growing high-order modes is re- 
stricted to the central 250 pc. Only the m=2-mode expands after 
t ~ 20 » 300 Myr to distances up to 1 kpc, but still with am- 
plitudes much smaller than in the center. In the nuclear region 
both, low- and high-order modes, reach a non-linear saturation 
level of about 5-20%. 



From a technical point of view, the conservation of vari- 
ous quantities is quite well fulfilled. At the end of the sim- 
ulation (after ~ 1.8 ■ 10 5 numerical steps) the mass is con- 
served within machine accuracy of double precision numbers 
(~ 10~ 16 ). Energy is conserved to better than 5 • 10~ 5 and an- 
gular momentum better than 2 • 10~ 6 (the numbers denote the 
relative deviations from their initial value). This simulation was 
performed without artificial viscosity. 



3. Results 

In this section we present the results for dusty disks. First, 
we discuss the properties of a reference model in Sect. 13.11 
Afterwards the influence of various parameters is shown in 
Sect. El 



3.1. The reference model 
3.1 .1 . The start model 

Our initial models are motivated by the nuclear region of Ml 00 
(NGC 4321). We adopted a total gas mass of 4.7 ■ 10 8 M© dis- 
tributed exponentially with a scale length of 300 pc within a 
radial range of R[ n = 30 pc to R out = 3 kpc. Thus, the central 
surface density of the gas is about 750 M0pc~ 2 . The small 
disk scale length mimicks a central concentration of (molecu- 
lar) gas observed in many galaxies. 

The rotation curve was selected as a combination of a cen- 
tral rigid rotation and an outer flat rotation curve according to 
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Fig. 4. Radial profile of the initial rotation curve and the differ- 
ent contributions to the gaseous azimuthal velocity of the ref- 
erence model: rotation curve, i.e. circular velocity v c (solid), 
azimuthal velocity of gas v rt >t (triangles), pressure contribu- 
tion to rotational equilibrium vp (dot-dashed), halo (and stellar) 
contribution to rotational equilibrium Vhaio.stars (short-dashed) 
and contribution form the self-gravity of the disk Vd; s k (long- 
dashed). 

Eq. dl 31 . The velocity at infinity, v ro , was set to 178 km s _1 . 
The transition parameter n, was selected to be 10, by this re- 
sulting in a fairly sharp transition at the radius =100 pc. 
This rotation curve corresponds to a total dynamical mass 
M d (R) ~ vl(R)R/G (including all components) of 7.3 • 1O 8 M 
within the central 100 pc. In the region of rigid rotation the rota- 
tion period is about 3.5 ■ 10 yrs. It increases outwards reaching 
1.7 ■ 10 7 yrs at the half-mass radius of the gaseous component 
at R * 500 pc. 

According to the chosen rotation curve, mass profile and 
equation of state (polytropic with y g = 5/3), the minimum 
value of the Toomre parameter is reached at a galactocen- 
tric distance of about 440 pc. The constant in the equation of 
state (for the gaseous phase) was selected to yield a minimum 
Toomre parameter of £>min = 1 .54. This choice corresponds to 
sound speeds between 4 and 1 1 km s _1 within the central kpc 
(the higher value is reached in the center). 

The different contributions to the azimuthal velocity of the 
gas as well as the rotation curve v c are shown in Fig.|4j The con- 
tribution vp of the pressure to the initial rotational equilibrium 
velocity is very small. It reaches a maximum of 8.5 km s _1 at a 
distance of 450 pc from the center. Compared to the azimuthal 
speed v rot of the gas vp is negligible as the almost identical 
values of the azimuthal speed and the rotation curve demon- 
strate. The radial force attributed to the gravitational potential 
is dominated by the contribution of the dark halo plus stellar 
disk/bulge. The self-gravity of our live disk components (gas 
& dust) gives rise to a maximum rotation velocity of about 50 
km s . This value still exceeds the pressure contribution by far, 
whereas it is small compared to the "halo" contribution made 
out of stars and dark matter. Hence, the (gaseous) disk is sub- 
maximal by one order of magnitude. 

Assuming that the central region is dominated by stars (i.e. 
the dark matter mass fraction is small there), the velocity dis- 
persion cr s of the stars must be a factor of 10 larger in or- 
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Fig. 7. Temporal evolution of the Fourier amplitudes of the m = 
0, 1,2,6,8, 10-modes of the dust component in the reference 
model. The time unit is 1.5 ■ 10 7 yrs. 



der to form a Toomre-stable maximum disk (their Toomre pa- 
rameter scales with Q ~ cr,/Z). Observations of the central 
stellar velocity dispersions report such large values of about 
o-j ~ 100- 150 km s~' (e.g. Heraudeau & Simien ll998> . Hence 
even a maximum stellar disk would be as Toomre-stable as the 
gaseous disk and more likely it would be dynamically hotter, 
i.e. more stable. On the other hand the scale height h of such 
a stellar disk is fairly large. Estimating h from h = ct 2 s /(7tG1,) 
yields for <x s = 100 km s -1 and Z = 1000 M pc~ 2 (the value 
of a maximum stellar disk at a galactocentric distance of about 
500 pc) a scale height of h ~ 750 pc. Hence, the stellar com- 
ponent is not flat at all within the central component. Together 
with the large Toomre parameter of the stellar component this 
justifies the treatment of the stellar "disk" as a "background" 
potential instead of a dynamically live component. 

In our reference model we set the dust mass to 2% of the 
gas mass. This choice is close to the upper value reported for 
the average dust-to-gas ratio in galaxies. On the other hand, the 
central regions of galaxies are metal-enriched compared to the 
mean galactic values. Taking the observed relation between the 
dust-to-gas ratio r and the metallicity into account (e.g. Issa et 
al. 119901 . it is reasonable to assume a larger r in the central 
regions than on average. 

Since we treat the dust as a pressureless component, the 
initial equilibrium azimuthal speed of the dust is equal to the 
circular speed v c , whereas the azimuthal velocity of the gas 
is slightly lower due to the pressure contribution to the radial 
forces. The velocity difference between both components, how- 
ever, is only about 0. 1-0.2 km s _1 within the central kpc. 

Though the choice of our initial model was guided by ob- 
servations of Ml 00, it should be kept in mind that the follow- 
ing simulations are not meant to yield a model of M100 2 , but 
to study a "typical" central galactic region. 

3.1 .2. Evolution of the reference model 

Fig. |5] shows images of the gaseous component at different 
times. During the first 300 Myr (t ~ 20 ~ 300 Myr) only very 



2 E.g. the weak bar is not included in our models. 
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Fig. 5. Surface density of the gas component of the reference model at different times. A white area corresponds to surface 
densities of 10 2 Mo pc~ 2 or lower, whereas black means 10 3 Mo pc~ 2 or higher. The length unit is kpc and the time is given in 
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Fig. 8. Temporal evolution of the Fourier amplitudes of the 
m = 0, 1, 2, 6, 8, 10-modes of the gas component in the refer- 
ence model. The time unit is 1.5 • 10 7 yrs. 



weak features are discernible in the gas phase. Lateron, a multi- 
armed and very patchy structure is formed which is basically 
concentrated to the central 400 pc. After 390 Myr a small re- 
gion close to the center becomes highly unstable. Since we did 
not include any processes like star formation or stellar feed- 
back the growth of this clump is not stopped. As a result the 
timestep in our simulations decreased to less than a year. At 
this value we stopped the numerical simulation. The formation 
of such clumps, however, is a generic behaviour of our simula- 
tions. When we varied the initial conditions by using different 
perturbations, clumps were still formed on a similar timescale 
(~ 500 Myr). 

In comparison to the gas, the dust becomes unstable at an 
earlier time (Fig. [5}: already after 100 Myr (f ~ 7) structures 
are visible in the central 100 pc. These structures are character- 
ized by a large contrast, i.e. thin, but (relatively) dense regions, 
separated by broad areas devoid of almost any dust. The po- 
sitions of the dust peaks are strongly correlated with those of 
the gas, but sometimes with a small offset in the positions. The 
absolute values of those peaks are not well correlated with the 
corresponding peak values of the gas density. This leads to a 
large scatter in the gas-to-dust mass ratios. Morphologically, 
the dust seems sometimes to be organized in rings instead of 
large-scale spirals like the gas. Some of these rings are linked 
by dust lanes or small arcs. 

A more quantitative description of the structure formation 
is given by the Fourier amplitudes. The dust component shows 
a fast exponential growth of all components with a slight dom- 
inance of the high-m modes (Fig. [7). Different to the single- 
component model discussed in Sect. I2.7l the dust modes begin 
to grow immediately at t = 0. Their growth rates are a factor 
of 2 larger than those found in the single-component model. 
Additionally, the evolution is accompanied by a strong radial 
mass redistribution of the dust. After 150 Myr a saturation level 
of about 10~ 2 is reached. 

Though the dust has only 2% of the mass of the gas, it 
destabilizes the gaseous disk strongly (Fig. |8j. When dust is 
present, the modes of the gaseous phase commence to increase 
already at about t ~ 5 ~ 75 Myr. In the single-component simu- 
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Fig. 9. Temporal evolution of the spatially resolved Fourier am- 
plitudes for the m=2-mode (upper diagram) and the m=8-mode 
(lower diagram) for gaseous component of the reference model. 
The amplitudes are calculated according to Eq. Jl 81 integrating 
over small radial annuli defined by the numerical grid resolu- 
tion. The time unit is 1.5 • 10 7 yrs. 



lation the gaseous disk has been stable until a time t ~ 10 ~ 150 
Myr (cf. Fig. 0. This difference becomes even more obvious 
when comparing the spatially resolved Fourier amplitudes of 
the simulations with and without dust (Figs.|3]and|9|i: the grow- 
ing region is for both simulations restricted to the central area, 
but in the dusty model the growth sets in much earlier and with 
a larger growth rate. 

Another difference to the single-component model con- 
cerns the dominant modes. Whereas for the purely gaseous 
model the high-m modes are dominant, all modes grow more 
or less at the same rate in the reference model. The superpo- 
sition of all these modes with nearly equal amplitudes results 
in the patchy structure visible in the spatial distribution of the 
components. 

The growth of the Fourier amplitudes of the gas shows four 
stages (Fig. [3}. In the first stage (until t ~ 5 ~ 75 Myr), the gas 
is almost unaffected by the presence of dust and remains on its 
initial perturbation level. In the second stage (until t ~ 10 » 
150 Myr), the gas reacts on the grown instabilities in the dust, 
i.e. the formed very dense dust lanes, by becoming unstable, 
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too. This phase of growth changes, when the dust reaches its 
saturation level. In the following third stage, the instabilities 
within the gas are still growing but on a longer timescale. The 
last stage is reached, when the gas phase saturates (which did 
not happen until the end of this simulation). 

The direct influence of the frictional force is visible when 
considering the velocities. E.g. the residual azimuthal velocity, 
which is normalized to the initial azimuthally-averaged equi- 
librium rotation speed, shows a small but systematic deceler- 
ation of the dust component. Because the dust is treated as a 
pressureless phase, its rotation speed exceeds that of the gas. 
Therefore, the frictional force leads to a deceleration and, by 
this, to an azimuthal velocity which is 0.05 kms -1 smaller 
than the circular speed. The dust gets out of rotational equilib- 
rium resulting in a steady dust inflow. This radial flow becomes 
larger with increasing coupling of the dust to the gas. Hence, 
the m = 0-Fourier amplitude (which characterizes the radial 
mass reconfiguration) varies stronger when increasing the cou- 
pling strength, while the higher Fourier amplitudes (m > 0) are 
less affected. 

The angular momentum removed from the dust is absorbed 
by the gas which is accelerated in azimuthal direction. Due to 
the large gas-to-dust ratio the acceleration is almost negligible 
and the radial structure of the gas is much less affected. The 
specific angular momentum of the dust component decreases 
by 6.6% during 450 Myr. The overall relative conservation 
of angular momentum is about 2 ■ 10~ 6 , i.e. the simulation is 
as accurate as the single-component model. The total energy 
decreased by 0.8% due to the dissipative nature of the fric- 
tional force. Though this decrease is small, it exceeds by far 
the intrinsic energy uncertainty of 5 ■ 10~ 5 found in the single- 
component simulation. Thus, the code is sufficiently accurate 
to deal with the implemented dissipation rates. Most of the en- 
ergy dissipation takes place when the clumps are formed and 
the system becomes highly non-linear. 

3.2. Parameter studies 

In a small set of parameter studies we investigated different as- 
pects of the dust's influence on nuclear gaseous disks. First, 
we compared in more detail the different dust-gas coupling 
schemes (Sect. l3~2~Tl . Then we varied the gas-to-dust mass ra- 
tio (Sect. l3~2~2l . In the sections B. 2. 3l and l3. 2.41 we investigate 
the influence of the Toomre parameter and the equation of state 
of the gaseous component. Finally, different minor aspects are 
summarized in Sect. 13.2.51 



3.2.1. Friction schemes 

In a first series of simulations we compared the different cou- 
pling schemes Eqs. © and (I12> applied to the reference model. 
Our simulations show that - independent of the scheme applied 
for the dust-gas coupling - the disks become more unstable 
when dust is present and the frictional timescale becomes not 
shorter than 10 4 yrs. 

In case of a constant frictional timescale according to Eq. 

we investigated four values of t^: 10 8 , 10 7 , 10 6 and 10 5 yrs. 



h — ' — i — 1 — r 



dust 




Fig. 10. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the dust component for different 
coupling timescales (r^ = 10 5 yrs (open boxes), = 10 6 yrs 
(plus), t d = 10 7 yrs (triangles) and 10 8 yrs (filled boxes)) and 
the reference model (thin disk friction; solid line). The dust-to- 
gas mass fraction is 2%. The time unit is 1.5 • 10 7 yrs. 



In case of = 10 8 yrs the frictional timescale is of the order 
of the rotational period at the outer boundary of the grid at 3 
kpc. This timescale exceeds the rotation period at the half-mass 
radius by a factor of 5.9 and the rotation period at the inner edge 
by a factor of 28.6. The dominantly growing modes of the dust 
component are multi-armed modes near m — 8. The growth 
of the two-armed mode is delayed. The dust modes grow very 
fast: already after 50 Myrs (f ~ 3.5) they reach their saturation 
level of about 1% (Fig. 1101. A weak coupling between gas and 
dust, i.e. a large t</, results in an almost decoupled dust disk 
which is extremely unstable due to its lack of pressure support 
on small scales. 

When Trf is reduced by a factor of 100, the frictional 
timescale is shorter than the dynamical timescale and gas and 
dust are strongly coupled. In that case the growth of pertur- 
bations in the dust component is reduced, but still the dust 
evolves to saturation within the first 10 8 years (Fig. I10i. The 
thin disk friction of the reference model grows on an even 
longer timescale, corresponding to a frictional timescale be- 
tween 10 5 and 10 6 yrs. The model with the shortest fixed fric- 
tion timescale of t,/ = 10 5 yrs differs qualitatively from the 
other simulations: whereas the other calculations exhibit an im- 
mediate growth of modes, the t</ = 10 5 -model remains stable 
for the first 120 Myr (f ~ 8) before it also develops instabil- 
ities. The coupling to the dynamically warm gas becomes so 
strong that the growth of instabilities in the dust component is 
efficiently suppressed. Independent of these differences, the fi- 
nal saturation levels for all friction strengths are quite similar 
of the order of 1%. 

In Fig. ^2 the growth of the dominant m = 8-modes of 
the gaseous phase are shown for different coupling strengths. 
Similar to the dust component, the evolution of the gaseous 
phase strongly depends on the frictional timescale. In compar- 
ison to the single-component model, friction strongly destabi- 
lizes the gaseous component, though its mass is a factor of 50 
larger than that of the dust. Also for the two more realistic treat- 
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Fig. 11. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the gas component for different 
coupling timescales (t^ = 10 6 yrs (plus) and 10 8 yrs (filled 
boxes)), the reference model (thin disk friction; solid line), the 
single-component model (no friction; dashed line) and a model 
with the "thick disk friction" (open boxes). The dust-to-gas 
mass fraction is 2%. The time unit is 1.5 • 10 7 yrs. 



ments of the dust friction, the thin and the thick disk limits, the 
evolution of the gas phase deviates significantly from that of the 
single-component model. In the thin disk limit, the gas starts 
to form structures 25% earlier than in the single-component 
model. If the frictional timescale becomes even much shorter 
like in the "thick disk" limit, the gaseous disk is even stabilized 
by the dust (only after t — 20 ~ 300 Myrs, the disk starts slowly 
to develop growing modes). In that case, the friction acts like a 
very high viscosity which slows down any growth of instabili- 
ties. 

The more physical approaches for the dust treatment are 
characterized by frictional timescales depending on the local 
gas and dust properties. In case of a thin dust disk, the frictional 
timescale is given by Eq. d!2i . i.e. the dynamical timescale. If 
we compare this ansatz with the simple approach of a constant 
frictional timescale, we find that the thin-disk scheme is very 
similar to a short frictional timescale between 10 5 and 10 6 yrs 
(Fig. llOl i. whereas the treatment according to Eq. gives an 
even shorter timescale. 

The coupling between the gas and dust component can 
be studied by comparing the evolution of the global modes 
of both components. If one applies Eq. the coupling be- 
tween gas and dust is very strong and the modes have nearly 
identical temporal evolution. If the frictional timescale is given 
by the dynamical time according to Eq. ill 2b the coupling is 
less strong (see Figs. 1101 and II Q . The Fourier amplitudes of 
the dust exceed then those of the gas throughout the calcula- 
tion. Especially in the phase of strongest growth, the dust is 
clearly ahead of the gas. This indicates that the instability of 
the gaseous disk is enhanced by the dust. 

3.2.2. Gas-to-dust ratio 

As a next step we varied the dust-to-gas mass ratio r = Md/M g 
from 0.5% to 20%. For all these models the m = 8-mode is 
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Fig. 12. Temporal evolution of the logarithmic Fourier ampli- 
tudes of the dominant m = 8-mode of the dust component for 
different dust-to-gas mass fractions: 0.5% (open boxes), 1% 
(crosses), 2% (reference model, solid line), 10% (plus), 20% 
(triangle) and the purely gaseous model (dashed line). The time 
unit is 1.5 ■ 10 7 yrs. 
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Fig. 14. Radial profile of the surface density of gas (solid) and 
dust (dashed) along a radial line at t — 5 for the 10% dust frac- 
tion model. The surface densities are given in 10 3 M0 pc~ 2 , the 
radius in kpc. 



dominant. A comparison of the corresponding Fourier ampli- 
tudes with those of the purely gaseous model of Sect. l2.7l shows 
that the critical dust-to-gas ratio r c for dust becoming dynam- 
ically unimportant is about 1% (Fig. I12>. For larger amounts 
of dust the destabilization of the gaseous phase becomes much 
stronger. For the reference model's value of r — 0.02 the insta- 
bility sets in 50 Myrs earlier than in dustless or low-dust models 
while the modes remain on their initial value for a latency pe- 
riod of 75 Myrs. Increasing r to 10% reduces the latency time 
to almost zero. The growth rates increase by a factor of 3-4 and 
the saturation level reached already after 30 Myrs is larger by 
at least one order of magnitude. 

Fig.ll3lshows the gas and dust distribution of a model with 
a large dust-to-gas ratio of 10% in its saturation stage. The 
weak structures visible in the reference model become more 
pronounced as a comparison of the gas distributions shows (cf. 
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Fig. 13. Spatial distribution of the surface density and its perturbations (normalized to their initial values) of the gas and dust 
component after t — 6 ~ 90 Myrfor a model with a 10% dust mass fraction ( relative to gas): E g ( upper left; white: 10 1 ' 5 Mq pc" 2 , 
black: 10 3 Mqpc~ 2 ), (upper right; white: 10 0,5 Mqpc~ 2 , black: 10 2 Mqp>c~ 2 ), AS,gj/Y, g j for gas (lower left) and dust (lower 
right). The latter are coded as follows: areas devoid of material are white, areas with a density enhancement of a factor of 2 or 
more are black. The grey area seen at the outer edges correspond to no deviation from the initial surface density. 



also last images in Figs.|5]and|6li: the patchy multi-armed struc- 
ture is more emphasized due to larger arm-interarm variations. 
Along the arms strong surface density variations exist. Some 
arms are interrupted by low density areas. Many spirals are not 
smoothly curved, but they show wiggles. Some arms seem to 
merge with others. The dust distribution is highly correlated 
with the gas distribution. However, the contrast between arm 
and interarm regions is larger for the dust than for the gas. The 
surface densities of the dust vary by about one order of magni- 
tude, whereas the contrast of the gas component is usually less 
than a factor of 2. (Fig. I14l i. The variations of the dust com- 
ponent along the arms are smaller than those of the gas. The 
structures formed in the dust are also thinner than those of the 
gas. Though there is a tight correlation between the positions 
of the maxima of the gas and dust phases, there is no correla- 
tion between their maximum amplitudes. The dust distribution 
is characterized by a more cellular appearance compared to the 
spiral-like morphology of the gas. Inspecting the density per- 
turbations stresses this point (see lower diagrams in Fig. ll3l >. 

Another interesting aspect is that the dust is often located 
at the boundaries around peaks of the gas distribution, prefer- 
entially at the inner boundary. Examples are the dust peaks at 



R = 125, 155 and 210 pc in Fig.^] Other dust peaks are in 
regions with no or only weak gaseous density enhancements 
like those at R = 150 and 170 pc. And, of course, there are dust 
peaks at the same locations as those of the gas, e.g. at R = 130 
and 190 pc. From that variety of locations it is clear that there 
is no simple straightforward correlation between the gas and 
dust mass distribution. Accordingly, one expects large spatial 
variations of the dust-to-gas ratio. 

3.2.3. Toomre parameter 

The Toomre parameter usually turns out to be a key indica- 
tor for the stability of disks. Therefore, we performed a series 
of simulations in which we studied dynamically cold, warm 
and hot disks, i.e. we varied Q from 1.1 to 3.0. The dust-to- 
gas ratio was set to the value of 2% identical to the reference 
model. As expected from single-component simulations Q has 
a strong impact on the structure formation of the gaseous phase 
(Fig.ll5>. Dynamically cold disks (Q close to unity) are highly 
unstable, whereas hotter disks become more and more stable. 
The saturation levels are much larger in case of cold disks (for 
which the non-linear regime is quickly entered). On the other 
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Fig. 17. Spatial distribution of the surface density perturbations ( normalized to the their initial values) at t — 16 ~ 240 Myr of 
the Q — 1 .2-disk: gas component (left), dust ( right). Areas devoid of material are white, areas with a density enhancement of a 
factor of 2 or more are black. The grey area seen at the outer edges correspond to no deviation from the initial surface density. 




Fig. 15. Temporal evolution of the logarithmic Fourier ampli- 
tudes of the m = 8-mode of the gas component for different 
initial (minimum) Toomre parameters of the disk: (2=1.1 (tri- 
angle), Q = 1.2 (x), Q = 1.54 (reference, solid line), Q = 1.8 
(dashed line), Q = 2.0 (boxes) and Q = 3.0 (dot-dashed). The 
time unit is 1.5 • 10 7 yrs. 




Fig. 16. Temporal evolution of the logarithmic Fourier ampli- 
tudes of the m = 8-mode of the dust component for different 
initial (minimum) Toomre parameters of the disk: Q — 1 . 1 (tri- 
angle), Q = 1.2 (x), Q = 1.54 (reference, solid line), Q = 1.8 
(dashed line), Q = 2.0 (boxes) and Q = 3.0 (dot-dashed). The 
time unit is 1.5 • 10 7 yrs. 



hand, even for large Q dusty gaseous disks are more unstable 
than single-component disks of the same Q stressing the desta- 
bilizing influence of the dust. 

Qualitatively, this can be explained by the evolution of the 
dust component whose growth rates are almost independent of 
the varied Q values (Fig. I16> . This is a consequence of the 
pressureless evolution of the dust which leaves the dust disk 
almost unaffected by the gaseous phase. There is just a small 
trend of the dust to become later unstable with increasing Q. 
Additionally, the loss rate of specific angular momentum of 



the dust component by transfer to the gas varies by a factor 
of 7.3 when Q is increased from 1.1 to 3.0. Both is caused 
by the larger pressure contribution to the rotational equilib- 
rium of the gaseous phase in case of dynamically hotter sys- 
tems. By this, the gas rotates more slowly, if Q is large, and 
the velocity difference between gas and dust becomes larger, 
too. For small pressure contributions to the azimuthal speed of 
the gas, the velocity difference between gas and dust scales like 
Av ~ ^ ~ K ~ Q 2 (K is the constant in the equation of state 
Eq. l|5}). As a result the frictional force is also enhanced ~ Q 2 



Theis Ch. and Orlova N.: Are galactic disks dynamically influenced by dust? 



15 




2 4 6 8 10 



time 

Fig. 18. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the gas and dust component for 
different equations of state (of the gas phase): y = 5/3, ref- 
erence model (gas: solid; dust: dashed) and isothermal (gas: 
crosses; dust: triangles). The time unit is 1.5 • 10 7 yrs. 
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Fig. 19. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the gas component for differ- 
ent initial perturbations: reference model (solid), another set 
of random numbers, otherwise identical to reference model 
(dashed), locally identical relative overdensities for dust and 
gas (boxes). The time unit is 1.5 ■ 10 7 yrs. 



and the onset of instability of the dust component is slightly 
delayed. 

It is interesting to compare the spatial matter distribution of 
dust and gas in more detail. The simulation starting with a cold 
disk of 2ini = 1.2 develops quickly a patchy multi-armed struc- 
ture in the gaseous phase. At the end of the simulation, when 
the dust reached already the saturation stage and the gas is in 
the non-linear regime of the growth of instabilities, the gaseous 
pattern is dominated by a multi-armed, irregular morphology 
outside 200 pc and a fairly regular, mildly changed distribution 
in the central part (Fig.^] left). In contrast to the gas, the dust 
is also rather irregularly distributed in the circumnuclear region 
(Fig -El right)- Outside this nuclear region, but still within the 
central 200-400 pc prominent multi-armed spirals are visible. 
In that region the dust component is strongly correlated with 
the gaseous phase. Moreover, in agreement with the reference 
model the structures in the dust are more pronounced than those 
seen in the gaseous phase. The difference in the locations of 
unstable regions between both components becomes more pro- 
nounced at earlier stages of the evolution, when the gas forms 
structures around a galactocentric distance of 400 pc, whereas 
the dust becomes mainly unstable inside the central 100 pc. 

Both can be explained in terms of the local Toomre param- 
eter Q g d of the components: Q g takes its minimum value of 1 .2 
at R ~ 400 pc growing to values of 1.5 at 200 pc. Therefore, 
the gas becomes preferentially unstable in a broad ring at about 
400 pc. On the other hand, there is no pressure support at all for 
the dust component. Hence, on small scales the dust is only pre- 
vented from collapse by the friction with the gas (and, by this, 
partly due to the gas pressure). On large scales, the dust (and 
the gas) can be stabilized by differential rotation. However, in 
the region of rigid rotation exists no differential rotation and, 
therefore, the dust can form quickly fragments or rings. 



3.2.4. Equation of state (gas component) 

In order to investigate the influence of the equation of state of 
the gaseous phase, we performed a simulation with an isother- 
mal equation of state (y = 1.0) instead of y = 5/3. The other 
parameters were kept identical to the reference model. Both 
simulations give nearly the same results as the growth rates of 
the dominant global amplitudes demonstrate (Fig. II 81 . Thus, 
the behaviour of the dusty disks does not strongly depend on 
the adopted equation of state, provided the minimum Toomre 
parameters are the same. 

3.2.5. Miscellaneous 

The influence of several "technical" parameters was studied by 
comparing the evolution of the dominant mode m = 8 with that 
of a reference model. 

Grid boundaries. As a first test we investigated the influ- 
ence of the radial extent of the grid by varying the inner and 
outer radial boundary by a factor of 0.5 and 2, respectively. 
The growth rates, however, remained basically unaffected by 
the variation of the boundaries. 

Grid resolution. In order to test the influence of the grid 
resolution, we performed a simulation with 686x686 grid cells. 
Both simulations show only minor differences, e.g. the evolu- 
tion of the dominant mode m — 8 is rather similar in both cases 
for dust and gas: the instabilities set in at about the same time, 
the growth rates differ slightly, i.e. they are marginally larger in 
the high-resolution model. The saturation levels, however, are 
identical. Therefore, we conclude that our results are indepen- 
dent of the grid resolution. 

Different seeds. The influence of the initial setup of per- 
turbations was tested in two ways: First, a different set of ran- 
dom numbers was used to determine the initial fluctuations. In 
a second simulation the initial perturbations of dust and gas 
were correlated by taking locally identical relative perturba- 
tions, i.e. an overdensity of gas corresponds to the same over- 
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Fig. 20. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the gas component for differ- 
ent artificial viscosity: reference model: no artificial viscosity 
(solid), C v i s = 1.0 (dashed), 2.0 (boxes), 3.0 (plus). The time 
unit is 1.5 ■ 10 7 yrs. 
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Fig. 21. Temporal evolution of the logarithmic Fourier ampli- 
tudes of the m = 8-mode of the dust component for different 
transition positions R^i of the rotation curve: reference model, 
^flat = 100 pc (solid) and Ra- dl = 200 pc (open boxes). The time 
unit is 1.5 ■ 10 7 yrs. 



density in dust (in the standard setup the perturbations were un- 
correlated). The global Fourier amplitudes were kept constant 
at a value of 10~ 6 . The dominant Fourier modes evolve very 
similar in all these simulations (Fig. 1191. The main difference 
is the end of the simulations, which was chosen to be the time 
when the numerical timestep drops below 1 yr. In all simula- 
tions this was initiated by the formation of a small scale gravita- 
tionally unstable clump. In all three simulations this happened 
at a time t ~ 30 ~ 450 Myr or later. Therefore, we conclude 
that the formation of the main structures is rather insensitive to 
the detailed setup of initial perturbations. 

Artificial viscosity. The use of artificial viscosity is a stan- 
dard tool in hydrodynamics in order to cope with shocks. 
Principally, shocks are expected in the late stage of the evolu- 
tion, especially during the non-linear or the saturation stage. On 
the other hand, artificial viscosity is formulated qualitatively 
similar to the terms describing friction, i.e. it couples to veloc- 
ity differences. Therefore, we did not apply artificial viscosity 
to the models shown so far in order to avoid a spurious mixing 
of effects induced by friction and artificial viscosity. 

However, we wanted to test the possible impact of (artifical) 
viscosity on our results. Therefore, we implemented the von 
Neumann-Richtmyer artificial viscosity described in Stone & 
Norman ( 1992 1 and modified the determination of the timestep 
in our code accordingly. The free constant C v i s (called C2 in 
Stone & Norman) describes roughly the number of cells over 
which a shock is smeared out. C v i s was varied from 1 to 3. The 
evolution of the global modes was basically unaffected by the 
application of artificial viscosity as the similar Fourier ampli- 
tudes in Fig. 1201 demonstrate. Hence, shocks seem to be less 
important for the structures formed in our simulations. This, 
however, might change when strongly barred nuclei are investi- 
gated as the simulations by Athanassoula ( 1992 1 have revealed. 

Transition to flat rotation. In the reference model the region 
which became unstable first was given by the transition from 
rigid rotation to the flat part of the rotation curve, i.e. at a dis- 
tance of Rs at = 100 pc (cf. also time t = 7 in Fig.|6j. Within that 



region the circular frequency is constant which means that a 
perturbation always acts on the same neighbouring gas packets 
and differential rotation can not pull apart those gas elements, 
by this preventing gravitational instability. Contrary to that, in 
the flat part of the rotation curve differential rotation stabilizes 
additionally the disk. 

Fig. [^compares the m = 8-modes for a rotation curve with 
a transition to the flat part at 7?g at = 200 pc instead of 100 pc as 
in the reference model. In agreement with the simple consid- 
erations of the previous paragraph the dominant mode grows 
faster: after a short initial period the region of instability ex- 
tends beyond 100 pc making a larger area unstable, by this en- 
hancing the growth rate of the m = 8-mode. Compared to other 
modes the dominance of the high-m modes is even more pro- 
nounced, because they affect a larger area and, thus, more mass. 
The saturation levels, however, are nearly identical in both sim- 
ulation which can be understood as a consequence of the only 
small mass fraction (< 10%) which resides in the area between 
100 and 200 pc. 

4. Discussion 

4.1. Preliminary remarks 

From previous works it is well known analytically as well 
as numerically that adding a cold component can destabilize 
galactic disks (e.g. Jog & Solomon 1984, Orlova et al. 2002 1. 
In our example the cold component is the dust phase, whereas 
the dynamically warm or hot component is the gas. A quick 
(but dirty) first estimate of the stability of the components can 
be done by inspecting the Toomre stability parameter, i.e. sta- 
bility is provided if Q exceeds unity: 

0~K 

Q = ^ > 1 • (20) 

cr is here the velocity dispersion of the component of interest, 
e.g. the sound velocity in case of a gaseous phase. This crite- 
rion is exact only for axisymmetric perturbations in flat, single- 
component gaseous disks. Qualitatively, however, it holds for a 



Theis Ch. and Orlova N.: Are galactic disks dynamically influenced by dust? 



17 



large variety of systems, if the RHS is replaced by a slightly 
larger numerical value. E.g. in case of stars one has a value of 
3.36/7T (BT87) or for non-axisymmetric modes one gets in case 
of flat rotation curves V3 (Polyachenko et al. U997> . 

For a two component system the Toomre criterion gives 
only a good estimate of the stability in some limiting cases in 
which the dynamics of the components can be well separated, 
i.e. the meaning of the dynamically active mass or surface den- 
sity S and its stabilizing pressure or velocity dispersion <x are 
well defined. This situation holds for some of the systems under 
investigation here, e.g. they can be described as cold dust disks 
embedded in a dynamically warm or hot gaseous disk which 
is in rotational equilibrium. The Toomre parameter is then de- 
termined by the external rotation curve, the dust distribution 
2 = Erf and its velocity dispersion. Since the dust is assumed 
to be pressureless, i.e. <x = the Toomre parameter vanishes. 
Thus, without other stabilizing processes the dust disk would 
become gravitationally unstable to all perturbations with short 
wavelengths 3 . 

In case of a weak coupling between gas and dust, the dust 
disk would evolve like a two-component disk. Especially the 
dust disk would become violently Jeans-unstable, provided the 
dust remains pressureless as assumed by Noh et al. (1991 1. 
The latter would be the case, if the increase of velocity dis- 
persion of the dust phase can be compensated by dissipative 
processes like inelastic collisions of dust particles. As a re- 
sult the dust might form dense regions, eventually dominated 
by dust. However, there is no evidence for the existence of 
such pure "dust balls" which means that there is either a suf- 
ficiently strong coupling between dust and gas or there is no 
efficient dust "cooling" (or both of them). In the case of proto- 
planetary disks Weidenschilling ( 1980) and Cuzzi et al. (1993) 
stressed the importance of turbulence which limits the den- 
sity in the dust layer, by this preventing the formation of dust- 
dominated regions. Similarly, small-scale turbulence in galactic 
disks might be an important (dynamical) heating mechanism, 
even in the absence of a drag force. On the other hand, in the 
limiting case of a very short frictional timescale the velocity 
of the dominant gaseous phase is simply imprinted on the dust 
component. The behaviour of the complete system is then again 
well described by that of a single component system character- 
ized by the gas parameters. 

The most interesting regime is that where the coupling be- 
tween gas and dust is intermediate and the dominant gaseous 
component is dynamically not too hot, i.e. not all instabilities 
of the gaseous component are completely suppressed. 

4.2. Is the dust component dynamically important? 

The simulations have revealed that the dust component is dy- 
namically important only, if the dust-to-gas mass ratio r ex- 
ceeds 1%. This value is close to values of r reported for the 
interstellar dust surrounding the solar system. In-situ measure- 
ments of the gas-to-dust ratio performed with the satellites 
Ulysses and Galileo in the heliosphere give a value of r ~ 1/94 

3 The same result can be derived from a more detailed stability anal- 
ysis of two-component systems (Jog & Solomon ll984t . 



(Frisch et al. 1999). If one corrects their sample (which was 
dominated by micron-sized grains) for smaller grains missed 
because of the interaction with the heliosphere r increases by a 
factor of 2 to ~ 2%, the value we used for our reference model. 

On the other hand, there are several dust-to-gas determina- 
tions from interstellar absorption lines towards different direc- 
tions (e.g. e CMa) which show a large scatter of r giving values 
down to only r ~ 0.2- 0.3% (Frisch et al. 119991 . Though the 
exact numbers strongly depend on details of the analysis (e.g. 
the mass distribution of the grains, the gas column densities, the 
adopted stellar abundances), the high local variability seems to 
be characteristic for the dust distribution. 

In external galaxies measurements of the dust-to-gas mass 
ratio are derived from IRAS observations. Young & Scoville 
(1991) give a mean value of r ~ 1/600. They explain the 
large discrepancy to the Milky Way (~ 1%) by the assump- 
tion that the bulk of dust might be cold (T < 30 K) radiating at 
wavelengths larger than 100 fim, whereas IRAS is mainly sen- 
sitive for "warmer" dust. Detailed determinations of the dust- 
to-gas mass ratio also suffer from difficult gas mass estimates. 
Especially, in the central regions where molecular hydrogen is 
often the dominant gas phase, the gaseous mass can only be 
inferred from CO measurements and a subsequent conversion 
to the H2 mass. Uncertainties in this conversion might easily 
accommodate a factor of two in the mass determinations. 

Though there are only a few dust-to-gas mass determina- 
tions, there is a clear close-to-linear correlation of r with the 
metallicity (e.g. Issa et al. [1990 Hirashita ll999> . In case of 
M5 1 the dust-to-gas ratio is about 2%, whereas for dwarf galax- 
ies like the Magellanic Clouds r is a factor of ten or more 
smaller. Edmunds & Eales ( 1998 ) estimated an upper total dust 
mass limit for spiral galaxies with a standard yield of p — 0.01 
to be 0.2% of the total baryonic mass fraction. Assuming a gas 
mass fraction of 10%, this gives an upper dust-to-gas ratio of 
about 2% averaged over the whole galaxy. In the central regions 
of galaxies, however, the effective yields are larger, e.g. Pagel 
( 1987 ) gives a value 4.5 times larger than in the disk, and the 
limits on r are accordingly higher. Therefore, dust-to-gas mass 
ratios of 1 % or slightly more seem not to be unusual high for 
normal spiral galaxies. 

Due to all these uncertainties, however, it is difficult to give 
a reliable global value for the dust-to-gas mass ratio in galactic 
centers. Local measurements yield dust-to-gas ratios exceeding 
the critical 1 % margin, whereas global ones favour smaller val- 
ues. However, both estimates might suffer from systematic ef- 
fects (and the differences of both are a hint for that). Therefore, 
it seems that dust-to-gas mass ratios of 1 % or more for normal 
galaxies cannot be ruled out from current data. Especially in 
the strongly obscured metal-enriched nuclear regions the dust 
mass fractions might exceed the locally determined values. 

It is also interesting to note that a strong place-to-place vari- 
ation of the dust-to-gas ratio is found in our simulations. This is 
especially remarkable, because we started with an almost uni- 
form dust-to-gas ratio (deviations of the order 10~ 6 ). Hence, 
even without an inhomogeneous dust production, large spatial 
gas-to-dust variations are produced. 
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Fig. 22. Temporal evolution of the azimuthally averaged sur- 
face density of the dust component at two neighbouring ra- 
dial cells (R ~ 140 pc) for the hot gaseous disk model with 
2min = 3.0. The time unit is 1.5 ■ 10 7 yrs, the surface densities 
are given in 10 3 Mq pc -2 . 



4.3. The friction instability 

An interesting difference between the single-component model 
and the dusty disks is that the dust forms ring-like perturbations 
in the central area, whereas the gas develops initially only high- 
m spiral modes (compare e.g. Figs. |2] and |6J. Obviously, the 
dust introduces an additional axisymmetric instability as the 
evolution of the azimuthally averaged surface densities of the 
dust component show (Fig.l22>. 

An axisymmetric instability due to a drag force has been 
predicted by Goodman & Pindor (2000). Their instability 
works without self-gravity, but needs the assumption that the 
drag force scales non-linearily (or more exact super-linear) 
with the surface density of the dust. Qualitatively, the mech- 
anism works as follows (see also Goodman & Pindor 2000 1: 
due to the friction there is a steady, but slow inflow of dust. If 
the frictional force scales super-linear with the friction be- 
comes stronger in areas with enhanced £,/ by this increasing the 
coupling to the gas and reducing the radial inflow. As a result 

starts to grow exponentially. 

As Goodman & Pindor (2000 1 pointed out, the drag de- 
scription of Noh et al. ( 1991 1 which we used here and which 
scales the frictional force strictly linear with the surface density 
does not allow for their instability mechanism. On the other 
hand, the self-gravity of the pressureless dust might play here 
a crucial role. The potential well of a ring-like perturbation 
acts like a super-linear frictional term which reduces the inflow 
velocity of the dust due to its gravity: dust flows from outer 
regions to the perturbed region, but cannot leave it inwards 
due to self-gravity. Additionally, the altered local gravity de- 
creases the density gradient of the gas component. Therefore, 
the pressure gradient is reduced and the difference of the ro- 
tation speeds of gas and dust becomes smaller. As a result the 
frictional force and the loss of angular momentum decreases, 
by this stabilizing the formed ring. That the dust's self-gravity 
is the driving agent is also supported by the fact that the rings 
form, when the spatially resolved Fourier amplitudes become 
non-linear. 
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Fig. 23. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the dust component for differ- 
ent treatments of the dust component: reference (pressureless 
self-gravitating dust; solid), no dust self-gravity (dashed) and 
non-vanishing pressure (boxes). The time unit is 1.5 ■ 10 7 yrs. 
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Fig. 24. Temporal evolution of the logarithmic Fourier ampli- 
tudes for the m = 8-mode of the gas component for different 
treatments of the dust component: reference (pressureless self- 
gravitating dust; solid), no dust self-gravity (dashed) and non- 
vanishing pressure (boxes). The time unit is 1.5 ■ 10 7 yrs. 

4.4. Variation of dust treatment 

As discussed at the beginning of this section, the implementa- 
tion of a persistently pressureless component, subject to self- 
gravity but without coupling to another dynamically warm 
component results in a very unstable disk. In order to inves- 
tigate the dependence of the coupled gas-dust system on the 
dust treatment, we performed two experiments, in which we 
modified the equation of motion of the dust by either describ- 
ing the dust with a small, but non-vanishing pressure term and 
by neglecting the self-gravity of the dust. 

These simulations demonstrate that the self-gravity of the 
dust is the main reason for the enhanced instability of dusty 
disks. In case of a disabled self-gravity of the dust, the dust 
"feels" only the more slowly developing potential wells of the 
gaseous phase. The Fourier amplitudes start to grow with a sub- 
stantial time delay, whereas they grow immediately in the ref- 
erence model (Fig. |23}. Once small perturbations in the gas 
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phase begin to form, the pressureless dust follows them imme- 
diately as the similar evolution of amplitudes shows (Figs. [23] 
and!24>. Due to the small mass fraction of the dust, the growth 
of the gaseous density perturbations is almost unaffected by the 
presence of a gravitation-free dust. In fact, the dust marginally 
delays the onset of growth of the gaseous perturbations due to 
the frictional force. 

Different to the self-gravity of the dust component, the 
treatment as a pressureless phase is less important, provided the 
dust is still dynamically cold. The growth rates of the Fourier 
amplitudes are merely identical, if an adiabatic equation of 
state is used for the dust as a model with Q m i n = 0.03 makes 
clear (Fig.l24>. The saturation levels do not vary with the dif- 
ferent treatments of the dust. 

4.5. Dynamics of dust grains 

For simplicity we investigated in this paper a dust component 
which is only coupled to the gas by a frictional force (and grav- 
ity). The collisional cross-section which enters the estimate of 
the frictional timescale in the thick disk limit (cf. Sect. 12.2.11 
and Eq. (|9}) is derived from the geometrical cross-section. 

Generally, it is assumed that dust particles are charged, and 
the sign of the charge depends among other parameters on the 
grain size. Such a charge leads to an additional Coulomb drag 
term for the interaction with ions. In case of fractional ioniza- 
tion (< 10~ 2 ) this contribution is small compared to the col- 
lisional term, whereas it can become a substantial contribu- 
tion in highly ionized areas (Draine 2003 1. In the thick-disk 
limit for calculating the frictional timescale A -1 an additional 
Coulomb term yields a reduction of the frictional timescale, 
whereas in the thin-disk limit A -1 is unaffected, because the 
frictional timescale is anyway dominated by the longer dynam- 
ical timescale. 

In addition to an enhanced Coulomb drag charged particles 
are also subject to a Lorentz force exerted by the interaction 
with the galactic magnetic field. For typical field strengths in 
the ISM of a few fiG Draine ( 2003 ) estimated a cyclotron pe- 
riod tb of the order of tb ~ 5 ■ 10 4 yrs. From a comparison of 
Tb with the collisional timescale t c of a few 10 3 yrs (cf. Eq. 
(II Q ), we find that the momentum exchange due to friction ex- 
ceeds the coupling to the magnetic field, provided we apply 
average ISM properties for estimating r c . On the other hand, in 
cool molecular clouds, the collisional timescale can be longer 
due to the reduced thermal velocities. Then the dust might be 
strongly coupled to the magnetic field. 

In addition to drag processes and a Lorentz force, the dust 
is also influenced by radiation pressure due to the (usually 
anisotropic) interstellar radiation field. Additionally, dust dy- 
namics is subject to recoil effects from photoelectric emission 
and photodesorption which can be of similar order as the ra- 
diation pressure. Draine (2003) estimated the drift velocities 
related to these radiation induced effects to be of the order of 
Vdrift ~ 0.04 kms -1 for a carbonaceous grain. Though this ve- 
locity is not very small compared to the relative velocities be- 
tween gas and dust in our reference model (~ 0.1 -0.2kms _1 ), 
we neglected radiation processes, because the related speeds 



are always small compared to the velocities reached during the 
growth of instabilities. Additionally, the determination of the 
interstellar radiation field would also require a consistent treat- 
ment of star formation and radiation transport in the ISM which 
is beyond the scope of this paper. The positional error resulting 
from the neglection of the drift velocities is only 4 pc within 
100 Myr. 

4.6. Miscellaneous aspects 

4.6.1 . Treatment of stellar component and dark matter 

In our simulations stars and dark matter were treated as a rigid 
component taken into account by its contribution to the rotation 
curve. By this, an exchange of angular momentum between the 
dusty gaseous disks and the stellar component or the halo is 
neglected. Such a coupling can be very important, especially 
for the large-scale evolution of galactic disks (e.g. Klypin et al. 
2002 1. On the other hand we focussed here mainly on the cen- 
tral kpc. In this region the mass of the dark matter is assumed 
to be small compared to the mass of the baryons (in case of 
normal galaxies). Thus, neglecting an interaction between the 
gaseous disk and a dynamically hot halo in the central area 
might be a less severe restriction. The situation is less clear for 
the stellar component. E.g. the perturbation exerted by a stellar 
bar will probably induce a strong m = 2-mode which might 
either dominate the appearance of the dusty disk or lead to a 
complicated mode-mode coupling between the dominant high- 
m-modes of the dusty/gaseous disk and the m = 2-bar mode. 

4.6.2. Formation and destruction of dust 

In this work we did not consider dust formation and destruc- 
tion processes. One reason was that we wanted to focus on the 
dynamical impact of the dust component on the stability of the 
gaseous disk and, thus, we wanted to isolate the influence as 
direct as possible. A second reason was that there is still no 
completely satisfying model for dust formation, though many 
aspects have been uncovered so far. From stellar outflows it 
is known that dust formation is strongly related to star forma- 
tion; other observations suggest that dust is formed in evolved 
stars or in supernovae events (e.g. Draine 2003 1. Similarly, the 
growth and destruction of dust shows a large variety of pro- 
cesses which are important for the evolution of the dust popu- 
lation. Therefore, the implementation of a detailed evolutionary 
model for the dust component seems to be out of reach at the 
moment. On the other hand, a simplified treatment of dust for- 
mation and destruction might be possible on the basis of a live 
stellar component which includes a simple description of star 
formation and stellar evolution. 

Due to the complexity of the processes governing the for- 
mation and destruction of dust, only rough estimates on the 
related timescales are possible. If we assume that dust forma- 
tion is mainly driven by stars (or supernovae), then the star for- 
mation timescale tsf provides a lower limit for the dust for- 
mation timescale rjf. For normal spiral galaxies the star for- 
mation rate is typically of the order of one or a few Mq /yr 
yielding a timescale of the order of (several) 10 9 yrs. A sim- 
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ilar lifetime was invoked by Draine & Salpeter (1979 1 in or- 
der to explain the large fraction of Si bound in dust grains 
(and not in the gaseous ISM). These lifetimes exceed other in- 
volved timescales like the revolution periods or the frictional 
timescales by far. Thus, dust formation and destruction is prob- 
ably less important for the overall dynamical processes pre- 
sented in this paper. On the hand, the secular evolution of the 
disk might be strongly affected, when the dust-to-gas ratio be- 
comes sufficiently large. Additionally, the spatial distribution 
of the dust-to-gas ratio might be influenced by locally concen- 
trated star forming events. 

These mentioned aspects (a live stellar disk, the influence 
of a stellar (mini-)bar and effects related to simplified dust for- 
mation and destruction model) will be part of a future paper. 

5. Summary 

We investigated the influence of a cold dust component on the 
evolution of galactic gaseous disks by means of 2-dimensional 
hydrodynamical simulations for flat disks. We focussed espe- 
cially on the question, whether a small contribution of dust is 
able to destabilize a gaseous disk and what kind of structures 
are formed. We coupled gas and dust by a drag force depend- 
ing on the relative velocities between both components and a 
frictional timescale. The latter is derived either only from col- 
lisions between gas and dust particles (thick disk limit) or from 
collisions followed by momentum equipartition in the gaseous 
disk (thin disk limit). Our initial models are composed of a sub- 
maximal exponential gaseous disk, a small admixture of dust 
and a rigid stellar and dark halo component. The rotation curve 
is assumed to rise linearily inside a radius of 100 pc becoming 
flat outside. 

From the evolution of the Fourier amplitudes we found that 
the higher-order modes are the dominant unstable modes in the 
purely gaseous simulation (in our case m = 8). Their growth is 
mainly restricted to the central kpc. An admixture of 2% dust 
(relative to the gas mass) destabilizes the gaseous disk. The 
dust component becomes non-linear after 100 Myr, followed 
by the gas 50 Myrs later. Different to the dust-free calculations, 
all modes have similar Fourier amplitudes. 

The formed structures of both, gas and dust, are rather 
patchy and multi-armed as expected from the similar Fourier 
amplitudes of the different modes. The dust component shows 
much more contrast between arm and inter-arm regions than 
the gas. The dust is spatially correlated with gas, but it does 
not exactly follow the gas. Peaks in the dust distribution are 
frequently at the inner edges of peaks in the gas distribution. 
Some dust peaks are also completely outside gas concentra- 
tions and some are exactly at the positions of gas peaks. This 
results in a large scatter of dust-to-gas ratios at different places. 
The dust develops also thin filaments which sometimes con- 
nect the arms. Therefore, the dust distribution has a more cel- 
lular appearance, whereas the gas develops a multi-armed spi- 
ral morphology. The dust is sometimes organized in a ring-like 
structure which is the result of an instability driven by the fric- 
tional force and the self-gravity of the dust (or adjacent gas). 

Below a dust-to-gas mass ratio r of 1% the dynamical in- 
fluence of the dust on the gaseous disk becomes negligible. 



This critical value is close to the observed mean value in nor- 
mal galaxies like the Milky Way. Since the dust-to-gas ratio 
scales linearily with metallicity, larger local values of r, espe- 
cially in the central galactic regions, seem to be reasonable. 
Such values are also in agreement with local gas-to-dust deter- 
minations. Since already a dust-to-gas ratio of 2% significantly 
affects the evolution of the disk, even the observed small dust 
admixtures are expected to have an impact on the dynamics of 
some galaxies (e.g. the dust-rich M51). For a 10% admixture 
of dust the gaseous component is completely destabilized. The 
growth rates are enhanced by a factor of 3-4 with almost no 
latency phase. The saturation levels reached after 30 Myr are 
substantially larger than in the low-r models. 

The Toomre parameter Q of the gaseous disk has almost no 
influence on the dust component, but it strongly controls the 
gaseous phase. As usual a larger Q gives more stable disks. 
However, the destabilizing influence of the cold dust was even 
found for a hot disk with Q — 3, though the saturation level is 
too small to be observable. 

The results are robust with respect to technical parameters 
like application of artificial viscosity, grid size and grid bound- 
aries. They are also independent on the adopted equation of 
state for the gas. For the dust treatment it is essential to take the 
self-gravity of the dust into account. The adopted equation of 
state of the dust is less important, provided the gas is treated as 
a dynamically cold component. 
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